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ABSTRACT 

A new statistical approach is presented to study the process of the thermal instability 
of optically thin unmagnetized plasma. In the frame of this approach the time evolu- 
tion of mass distribution function over temperature <p(T) is calculated. Function (p(T) 
characterizes the statistical properties of the multiphase medium of arbitrarily spaced 
three-dimensional structure of arbitrary (small or large) temperature perturbations. 
We construct our theory under the isobarical condition (P = const over space), which 
is satisfied in the short wavelength limit of the perturbations. The developed theory 
is illustrated in the case of thermal instability of a slowly expanding interstellar cloud 
(smooth scenario). Numerical solutions of equations of the statistical theory are con- 
structed and compared with hydrodynamical solutions. The results of both approaches 
are identical in the short wavelength range when the isobarity condition is satisfied. 
Also the limits of applicability of the statistical theory are estimated. The possible 
evolution of initial spectrum of perturbations is discussed. The proposed theory and 
numerical models can be relevant to the formation of the two-phases medium in the 
~ lpc region around quasars. Then small warm (T ~ 10 4 K) clouds are formed as the 
result of thermal instability in an expanded gas fragment, which is a product of either 
the star-star or star-accretion disc collision. 

Key words: hydrodynamics — instabilities — ISM: general — plasmas. 
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1 INTRODUCTION 

The interstellar and intergalactic optically thin plasma may 
be in a variety of thermal phases, depending on local heating 
and cooling processes and past history. The equilibrium of 
the thermal phases in the case of thermal and ionization bal- 
ance and in the case of constant pressure of the plasma have 
been widely discussed (for references see, e.g., Lepp et al. 
1985; McKee & Begelman 1990). More recent investigations 
are concentrated on the process of two-phases medium for- 
mation (Aranson, Meerson & Sasorov 1993) and dynamics 
of this medium (Elphick, Regev & Shaviv, 1992; Aharonson, 
Regev & Shaviv 1994). 

We focus our attention on the dynamics of the thermal 
phases generation and separation of mass over temperatures 
due to the the thermal instability process. We are interested 
in the evolution of the most rapidly growing modes which 
have time scales (and corresponding spacial scales) compa- 
rable with the characteristic cooling time t c . The growth of 
small amplitude perturbations of temperature and of den- 
sity of the optically thin medium under the action of exter- 
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nal heating and radiative cooling was studied by Weymann 
(1960) and Field (1965) in the linear approximation. The 
nonlinear regime is traditionally investigated by solving the 
exact or reduced hydrodynamical equations for the density, 
velocity and temperature perturbations of gas in space (see, 
e.g., Aranson et al. 1993). But, the detailed evolution models 
in these approaches can be only constructed by assuming the 
one-dimensional (ID) distributions. The attempts to solve 
this problem in two- and three-dimensions meet significant 
technical difficulties. 

In the present paper we propose a new statistical ap- 
proach to discribe the dynamics of the thermal instabil- 
ity and construct the theory for any (linear or nonlinear) 
stages in isobaric approximation. In the frame of this ap- 
proach we investigate the time evolution of the mass distri- 
bution function over temperature ip(T) [see equation (3.3)], 
which characterizes the statistical properties of the mul- 
tiphase medium. Such a medium may contain arbitrarily 
spaced three-dimensional perturbations of short scales. We 
assume the constant gas pressure (P = const) over volume 
even large contrasts of density p and temperature T exist. 
But, at the same time, the pressure is a time dependent 
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value. So, at any moment t the density and temperature of 
matter meet the relation p oc P/T. 

We show in Section 3 that in the case of the isobaric 
and nonuniform medium the temperature growth rate T of 
the matter with temperature T [see equation (3.12)] is 



T oc 



1 1 — T 

-c + (i--)cL 

7 7 T 



(1.1) 



Here C = C(p,T) = C(P/T,T) is the local cooling-heating 
rate of the gas, C and T are the mass average values of C and 
T, respectively [see equation (3.11)], and 7 is the gas adia- 
batic index. In the case of the uniform medium the equation 
(1.1) reduces to the form T oc —C One can see that the 
growth rate (1.1) depends on the local T only. Other quan- 
tities are common for any parts of the medium. They are 
the average over mass parameters, but are not independent 
variables. In this case T = T(T) the time evolution equation 
of the function <p(T) [see equation (3.15)] is of the Liouville 
type: 



0. 



(1.2) 



These coupled equations (1.1) and (1.2) determine the evo- 
lution of the medium on any (linear or nonlinear) stages. 

We illustrate our statistical approach in the case of 
a slowly expanding interstellar cloud when the medium 
smoothly evolves from the thermally stable stage to the un- 
stable stage. This evolution occurs through the critical point 
of marginal stability. We keep in mind the evolution of the 
initially 'warm ' (T ~ 10 4 K) gas because it can be in a sta- 
ble equilibrium in a very wide range of the external heating 
and ionization parameters. The initial stages of the thermal 
instability follow the linear theory of Weymann (1960) and 
Field (1965). On the subsequent nonlinear stages, when the 
linear theory is not applicable, we compare the numerical 
solutions of the equations (1.1) and (1.2), and ID hydro- 
dynamical solutions. The results of both statistical and hy- 
drodynamical approaches are identical in the isobaric limit, 
when the pressure is almost constant over the medium. We 
also estimate the limits of applicability of the statistical ap- 
proach. 

The hydrodynamical approach is discussed in Section 2. 
In Section 3 the statistical theory of thermal instability is 
constructed in the isobaric approximation. The smooth sce- 
nario of the thermal instability process is considered in Sec- 
tion 4. In Section 5 we discuss the numerical solutions of 
the equations of the statictical approach and hydrodynami- 
cal equations, and compare the results of both approaches. 
In Section 6 we discuss the obtained results and the limits 
of the constructed theory. In Appendix A we describe the 
numerical procedure of solution of the statistical equations. 
In Appendix B we present the results of the linear theory in 
the considered case of the smooth scenario. 



2 HYDRODYNAMICAL APPROACH 

The dynamics of a thermally unstable plasma can be de- 
scribed by the hydrodynamical equations for ideal gas with 
the additional cooling-heating term in the energy equation. 

dp 



dt 



+ pVv = 0, 



(2.1) 



dv „ 

P— + VP = 0, 



dt 



dP 



7" 



1 dt 



+ ■ 



-PVv + pC{p,T) = 0, 



P = *pT. 



(2.2) 
(2.3) 
(2.4) 



Here p, P, T, v are the density, the pressure, the tempera- 
ture, and the velocity of gas, respectively, R is the gas con- 
stant, p is the mean molecular weight, and 7 is the adiabatic 
index. For simplicity, we do not take into account the in- 
fluence of the ionization processes on thermal balance, and 
assume 7 and p to be constant. We assume fully ionized 
hydrogen plasma, setting 7 = 5/3 and p — 1/2 for numeri- 
cal estimations. The cooling-heating rate C is defined as the 
rate of energy losses minus the heating rate per unit mass. 
We assume that C depends on the local density p and tem- 
perature T of gas. In the simplest case of low density and 
optically thin plasma (Spitzer 1962) this rate is 



C(p,T) = P A(T) - H, 



(2.5) 



where A(T) and H are cooling and heating rates, respec- 
tively. The heating rate of the plasma can be determined 
by different mechanisms (e.g., photoabsorbtion, Compton- 
effect, heating due to cosmic rays and nuclear decays). For 
generality, we do not specify here these mechanisms and use 
a constant value for the heating rate in our models. In the 
simplest case of a stationary medium of temperature To and 
density po the criterion of marginal stablibility takes the 
form (Weymann 1960, Field 1965) 



d_ 

dT 



(?) 



= 0, po 



A(T ) 



(2.6) 



In the stable region the derivative is positive, in the un- 
stable region it is negative. The cooling function A(T) was 
calculated in a number of papers (e.g., Cox & Tucker 1969; 
Raymond, Cox & Smith 1976) for the different conditions 
in the interstellar plasma. For typical conditions of thermal 
interstellar plasma the marginal temperature To is about 
1Q A K. The function A(T)/T has a maximum at T = T . 
At T < To the plasma is thermally stable, and at T > To 
the plasma is unstable. In the presence of external ionizing 
(UV) radiation the function A(T) changes significantly (e.g., 
Krolik, McKee & Tarter 1981). 

In equation (2.3) we omit the thermal conductivity term 
V(ftVT), where k is the thermal conductivity coefficient. 
The thermal conductivity process stabilizes the thermal in- 
stability at the shortest wavelengths A & 2ir(hip/nRp) 1 / 2 , 
where n is the growth rate of instability (see the discussion 
in Section 6) . The modes of larger wavelength can be devel- 
oped till the nonlinear stages. For these modes the thermal 
conductivity process is important at the nonlinear stage, 
when the medium becomes very nonuniform and is divided 
into cold and hot phases. The thermal conductivity effec- 
tively stabilizes the further temperature growth of the hot 
phase and can produce the evaporation of the cold phase. In 
a subsequent paper, we will analyze the influence of thermal 
conductivity at the nonlinear stage of the thermal instability 
process. 
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3 STATISTICAL APPROACH 

We consider the cloud of plasma of mass M and of volume V . 
The equations (2.1)-(2.4) give us a precise description of the 
thermal and dynamical evolution of this matter. However, 
the solution of the problem in general case meets signifi- 
cant technical difficulties. We simplify the analysis propos- 
ing a new statistical approach which is constructed under 
isobaric approximation. The isobaric condition corresponds 
to a short wavelength limit 



A ^ Ap = 2irc s t c 



(3.1) 



when the pressure balance time ~ A/c s is shorter than the 
characteristic thermal time 

to = f. (3.2) 

Here c s = ("(P/p) 1 ^ 2 is the adiabatic sonic speed. The hydro- 
dynamical calculations confirm (see Section 5) that in the 
wavelength range (3.1) the isobaric approximation is satis- 
fied. The calculations also set the limits of applicability of 
the statistical approach. 

We build our theory using the temperature-dependent 
distribution of matter instead of the space-dependent one. 
We sum the masses of equal temperature and compose 
the distribution function <p(T) of mass over temperature. 
Namely, the mass of gas Am having the temperature in the 
range from T to T + AT is 



Third, taking the average in (3.9) we set T = PVp/MR 
[see equation (3.5)] instead of T in the right side and cor- 
respondingly (PV)p/MR instead of T in the left side, and 
find for the pressure derivative 



P = - (7 -l)fZ-4p. 



(3.10) 



Here and below the mass average of any function A(T) is 
defined as 



A = ip(T)A(T)dT. 



(3.11) 



Finally, substituting P from (3.10) into (3.9) we have the 
temperature growth rate of the matter with temperature T: 



7~ 1 M 
7 R 



(7 



—T 
1 £= +£ 

T 



(7- 



V 

DyT. 



(3-12) 



This equation follows from the mass (2.1) and energy (2.3) 
conservation equations. Instead of the Eulear equation (2.2) 
we use the isobaric condition P = C — const. Hydrody- 
namical calculations give us that P = C + 0(A 2 ) (see Ap- 
pendix B), but we ignore here the additional term 0(A 2 ) in 
the A ^ Ap limit. 

Note, that small velocities on small scales provide a 
scale-independent value of the velocity divergency, Vv ~ 
v\/\. One can substitute (3.10) and (3.12) into (3.8) and 
obtain explicitly 



Am = Mip(T)AT, (3.3) 
where the distribution function ip(T) is normalized to unity 



ip(T)dT = 1. 



(3.4) 



Now we show that the temperature growth rate T de- 
pends only on local T and the global (average) parameters of 
the medium. First we prove that for isobaric perturbations 
the pressure of the gas 



a V 

is determined by the mass average temperature 



T<p(T)dT, 



(3.5) 



(3.6) 



and the average density of the cloud (p) = M/V . Indeed, 
using equations (2.4) and (3.3) we can write 



f dm = r Mip(T)dT = R M [°° Tip{T)dT 

Jm P Jo P P Jo P 

In the isobaric case P = const the equation (3.5) follows 
from (3.7) directly. 

Second, from equations (2.1) and (2.4) we have for the ve- 
locity divergency 



Vv 



P T 
P + T* 



(3.8) 



Substituting Vv into equation (2.3) we obtain the relation 
between temperature and pressure growth rates 



(3.9) 



7 r It tj v 



(3.13) 



which is determined only by the local temperature T. 

In the isobaric case, according to (2.4) and (3.5) we 
have for density p = {p)T/T. It has to be set in the function 
£(p,T) everywhere. For example, the rate C, determined by 
formula (2.5), and its average value C, both of which we 
need to put into equations (3.12) and (3.13), have forms 



C = { P )(p}T-H. 



(3.14) 



After that we have all of the parameters (£, C, T) to substi- 
tute to formula (3.12) for the rate T. Also we have formulae 
for the pressure (3.5) and the velocity divergency (3.13) of 
the medium. 

The time evolution of the mass distribution tp(T) we 
describe by the Liouville type equation 



^ + ^(T^) = 0, 



(3.15) 



where the temperature growth rate T is given by formula 
(3.12). The equation (3.15) is nonlinear because the term T 
depends on the function ip(T) through the average values 
T and C. If we specify the initial conditions at the moment 
ti, the cooling-heating function C(p,T) and the law of ex- 
pansion of the medium V(t), then the equations (3.12) and 
(3.15) describe the thermal evolution of medium for t > U. 
The initial conditions are the distribution function tp(T, ti) 
and mean density {p)i. To solve the equations (3.12) and 
(3.15) we use a numerical technique (see Appendix A) in 
the general case of finite amplitude perturbations. 

The evolution equation for the relative temperature per- 
turbation 8T/T, where ST — T—T, can be constructed using 
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equation (3.12) for f and its average version for T. One can 



d_ /ST_ 
dt\T 



1 UzL-L 



7 RT \ T 



Under linear approximation we can write C(T) = C(T) + 
(dC/dT) p ST and find Z = £(T), because of ST = 0. In this 
case the equation for the relative temperature perturbation 
becomes 



d_ 

dt 



ST 
T 



7 RT V \dTjp J T ' 



(3.16) 



Here the values of C and (dC/dT) p are taken at T = T. 
For example, in case of the static uniform medium at the 
thermal equilibrium £(T) = we find the growth rate 



7" 



7 -RT 



(P) 



r 



\&tJ , 



(3.17) 



which coincides with the growth rate found by Weymann 
(1960) and Field (1965). For the cooling-heating function 
(2.5) it is reduced to 



7" 



tc 



TdA\ 

A dry 



(3.17') 



where t c is given by (3.2). From equation (3.17') one can 
obtain the criterion (2.6). 

Finally we note the following: 
1. One can construct uniquely the distribution function <p(T) 
from any T(r) summing the correspondent masses of vol- 
umes AVi with temperatures in range from T to T + AT, 



(p(T)AT = 



M 



For instance, consider a one dimensional case and the cosine 
profile of temperature, T(x) — To — STo cos (2irx/X). Then 
we can obtain inverse dependence x(T) for the interval < 
x < A/2 of monotonicity of T(x) and find the distribution 
function 



= 2Td*(T) 
^ v ' XT dT 



IT 



To 



(T - To) 2 
STo <T<To + ST , 



(3.18) 



where the mass average T 
dispersion a\ = (T 



T) 2 /T 2 



i(5T /T ) 2 at ST < T,. 



The function (3.18) is independent of the space period A 
and has two integrable infinities at T m in ~ To — STo and 
Tmax = T) + STo as shown in Fig. 1. 

2. The reconstruction of the space distribution T(r), p(r) 
and velocity field v(r) from a given <p(T) is an incorrect 
procedure. The recovery can be uniquely done in the one di- 
mensional case assuming additional restrictions on the tem- 
perature profile T(x). The temperature must be a periodical, 
monotonical for one half period (0 < x < A/2) and odd sym- 
metrical (relative x = 0) function. Keeping this restrictions 
one can find the implicit dependence T(x): 



A T., 



T(x) 



T'(p(T')dT'. 
Then for the density distribution we have 



p(x) = (p)T/T(x). 



(3.19) 



(3.20) 




-0.5 0.5 

(T-T )/ST 

Figure 1. The distribution function ip(T) in the case of the 
cosine profile T(x) = To + STo cos(2nx/ X) for different values of 
relative amplitude ST /T = 0.01, 0.5, 0.99. 



The velocity profile can be obtained by integrating equation 
(3.13): 

V 

v(x) — v(x) + —x — 



TO) /-£ 



V 



=T'-£(T')j <p{T')dT' + -x. (3.21) 



7 — 1 /i A 
7 RT 2 J 

The amplitude of the perturbed part v(x) of v(x) is propor- 
tional to the space period A, and 5(0) = v(X/2) — 0. 



4 SMOOTH SCENARIO 

We illustrate our statistical theory of thermal instability in 
the case of inertialy expanding uniform plasma cloud. The 
density of the cloud decreases with time, say, 

P u(t) cc i (4.1) 

Initially the plasma is assumed to be quite dense and cold. 
The temperature T u of the uniform plasma is less than the 
marginal value To , and the plasma is thermally and mechan- 
ically stable [see condition (2.6) and subsequent discussion]. 
In the case of a constant heating rate H the temperature T u 
increases up to the marginal temperature To. It happens at 
some moment to- At time t > to the plasma is unstable. 

Numerical models show (see Section 5.1) that the frag- 
mentation of cloud is only possible in the case of a 'slowly' 
expanding cloud for which the characteristic thermal time 
t c [equation (3.2)] is much shorter than expanding time to. 
In this case t c to the thermal (quasi) equilibrium of uni- 
form matter is established, C = 0. Then the evolution of 
temperature T u (t) follows from the equation 



A(T„) 



H_ 

pu 



A(To)(l) 3 . (4.2) 
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Near the marginal point To (where A 
find 

St, 



TdA/dT) one can 



(4.2') 



T u (t) = T (l + 3- 
to 

where St = t — to <C to. 

In the thermally unstable plasma the fastest growth rate 
has a condensation mode which corresponds to the condi- 
tions A £ Xp and P « const. In the linear regime the am- 
plitude A of this mode follows the equation A = nA At 
St <C to the growth rate n is small and linearly proportional 
to St: 



n(t) = 3b 



7- 



1 St 



(4.3) 



t c to 

where b = -T§ A"(T )/A(T ) is the positive quantity. The 
formula (4.3) follows from equation (3.17'), where the term 
in the brackets can be linearly approximated 1 — TA'/A = 
-T (T - T )A"/A at the limit T — To < T . In addition 
the time dependency of temperature (4.2') is used. Here and 
below we define the parameter t c — c 2 (To)/H at the moment 
t = to- On the linear stage ( Weymann 1960, Field 1965, see 
also Appendix B) at < St < r, the amplitude increases as 



A(t) = A expN(t), N(t) 



ndt - 



2 t, 



ISt 2 



t 



(4.4) 



and we must study the progressively inhomogeneous 
medium. Equation (4.1) is then correct for the average den- 
sity (p) only. The linear regime is limited by the time interval 
< St < t, where the end time 



2t t e ln(l/A]) 
36(7 - 1) 



(4.5) 



At St = t the integral N = ln(l/^4o) and the amplitude 
A — 1. At t > to + r the perturbations are large and 
grow nonlinear ly. In case of the small parameter t c /to <C 1, 
the time interval of the linear stage is also small: r/to ~ 

yJln(l/A )t c /t < 1. 

The scenario of the expanding cloud is able to follow 
the natural smooth evolution of the medium from the stable 
(t < to, T < To) to the unstable (t > to, T > T ) stages 
through the critical point (T = To). Other scenarii of a 
smooth transition through the critical point T = To can 
also be constructed. For example, we can take the stationary 
medium (p u — po — const) heated at a slowly increasing 
rate and obtain similar results assuming t c <C to and setting 
H(t) = Hoit/tof, Ho = poA(To). 



5 NUMERICAL RESULTS 

We use the expanding cloud scenario to numerically study 
the thermal instability process in the framework of statisti- 
cal isobaric theory. We find numerically the time evolution of 
the distribution function <p(T, t) and the evolution of corre- 
sponding characteristic parameters of this distribution. We 
also find the final mass fraction of cold and hot phases as a 
function of the parameters of the problem. To test the re- 
sults of the statistical approach we compare them with the 
results of hydrodynamical calculations in the following way. 
We use the obtained <p(T, t) to construct the T(x), p(x) and 
v(x) profiles. To do this we use the reconstruction procedure, 



which is discussed at the end of Section 3 [see equations 
(3.19), (3.20) and (3.21)]. After this we compare the con- 
structed profiles with the profiles resulting from the ID hy- 
drodynamical calculations under the identical parameters of 
the problem and the similar initial conditions. In all models 
we fix the heating rate H = const and set 7 = 5/3, p = 1/2. 
The expansion law (4.1) for the average density (p) of the 
plasma is assumed in the calculations. 

We use the model curve of the cooling function A(T), 
which mimics the base properties of the cooling processes in 
the interstellar plasma, 



where the function 8(T) is 



(5.1) 



9{T) = 



exp[d(l/c-r /T)], 



if T/To > c, 
otherwise. 



Here we assume c = 0.9 and d — 10.0. The exponential drop 
of 6(T) in the low temperature region T < cTo mimics the 
strong fall of A(T). The cooling curve (5.1) has the criti- 
cal point T = To = 10 4 K, where the criterion of marginal 
stability (2.6) is fulfilled. The instability growth rate coeffi- 
cient [see equation (4.3)] b = 1 in this case. We note, that 
the cooling curve (5.1) does not allow a stable equilibrium 
for a two-phases medium. 

We start our calculations at moment t = to. The initial 
relative perturbations of temperature A = (T max — T m m)/T 
and the expansion parameter to/t c = toH/c 2 (To) are free 
parameters of the models. The hydrodynamical models, in 
addition, are characterized by the initial space period Xo/Xp 
of perturbations. Here Xp = 2wc^(To) / H = 2nc s t c . 

The cosine profile of tempera- 

ture T(x) — To (l — y cos(27ra;/Ao)) is used as the initial 
perturbations in the hydrodynamical models. In the statis- 
tical models the corresponding initial distribution function 
<fi(T) is given by equation (3.18). The initial temperature 
dispersion erf (to) = A 2 /8 = A%/2. 



5.1 Statistical isobaric models 

The numerical method used to solve the equations (3.12) 
and (3.15) of statistical approach is briefly described in Ap- 
pendix A. 

We discuss, as a characteristic example, the detailed 
evolution of the S-model, with the initial parameters A = 
10 -4 and to/tc = 10 4 . The temperature dispersion a 2 -. — 
(T — T) 2 /T 2 increases with time at the t > to stage. The 
evolution of ar(t) is shown in Fig. 2. Initially or (to) = 
A/V8. The solid line corresponds to the nonlinear cal- 
culation, and the dashed curve was calculated under the 
linear approximation [equation (3.16)]. As it was shown 
in Section 4 [equation (4.4)] at the linear stage <xt(£) ~ 
c"T(to) exp [(t — to) 2 /tot c ]- A significant deviation is ob- 
served between the linear and nonlinear stages when the 
value of ot is large, <tt <) 0.1. The instability results in 
the formation of the two-phases medium at the moment 
t ~ to + T, where the time interval r is determined by equa- 
tion (4.5). 

In Fig. 3 (left- and right-hand panels) we show the tra- 
jectories of mass particles in (T,t)-plane. Each line is the 
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Figure 3. The evolution of temperatures of the fixed mass coordinates (solid lines) and average temperature T (dashed line) for the 
S-model. The mass coordinates of neighbouring solid lines differ by AM/M = 0.05. The upper and lower lines correspond to maximum 
and minimum temperature of plasma. The long time evolution is shown in the left-hand panel. Detailed evolution at right represents the 
nonlinear stage shown inside the box in the left-hand panel. 




t/t 

Figure 2. The evolution of the relative perturbation of temper- 
ature <tt for the statistical isobaric S-model. The solid line shows 
the evolution coming from nonlinear theory, and the dashed line 
shows the result in the linear approximation. 

trajectory of a fixed mass coordinate, and for each two neigh- 
bour lines li and £2 the following condition is satisfied: 

AM f T{t2) 

for any time t. The upper and lower lines correspond to max- 



imum and minimum temperatures of plasma. The dashed 
line shows the evolution of average temperature T(t). Note, 
that the time dependence of T(t) does not coincide with 
the evolution of temperature of the uniform medium fol- 
lowed from equation (4.2). In the linear stage of instability 
the lines are spaced closely. In the nonlinear stage, which is 
approximately shown inside the box in the left-hand panel 
and in more detail in the right-hand panel of Fig. 3, the 
lines disperse rapidly. One group of lines goes up to the high 
temperatures, and another group goes down to the low tem- 
peratures. A small part of gas of intermediate temperatures 
shows nonmonotonic behaviour in its temperature evolution. 

At the following time t > to + r the two-phase medium 
with hot (Th, ph) and cold (T c , p c ) phases is formed. A 
sharp boundary separates phases. Numerical calculations 
show that the mass fraction of the cold and hot phases tend 
asymptotically to final values. In the case of the S-model the 
final fraction of the cold mass is r\ — 0.86. The remainder 
1 — 77 = 0.14 is the final mass fraction of the hot matter. The 
temperature dispersion of this medium is 

4=(|-l) 2 (l-,)+(l-|)%. 

In the limit T c <C Th one obtains a simple relation a% = 
?7/(l — 77). For the case r\ — 0.86, we find <tt = 2.5, which 
very well agrees with our numerical result (see Fig. 2). 

Consider now the evolution of the parameters of hot 
and cold phases at the stage t — to <C to, when the expansion 
can be neglected. Assuming Th 3> T c one can substitute in 
equation (2.3) the density of hot gas ph ~ (p)o(l — 17), where 
(p)o = H/Ao, and obtain the equation: 

r fc = ( 7 -i)£jr(^(.j-i) + iV (5.2) 
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Figure 4. The final mass fraction of cold phase 77 of the statistical 
isobaric models as a function of the expansion parameter to/tc 
for the initial amplitudes A = 10 -2 , 10~ 3 and 10~ 4 (deltas, 
full circles and open circles). A pentagon and square show the 
final mass fractions of the hydrodynamical HI and H2-models, 
respectively, which have the initial amplitude A = 10 — 4 . The 
Hl-model and isobaric S-model (at to/t c = 10 4 and A = 10 — 4 ) 
have almost equal final mass fractions. 



The ratio A(T^)/Ao = 2 [see formula (5.1)] and we obtain 
Th ~ (2?7 — 1). The temperature of the hot phase increases 
to infinity and a two-phase medium can be formed in the 
case when more than half of the matter falls down to the 
cold phase, 77 > 1/2. Finally, at t > to + r, we have 

T fc ~(7-l)ii#(277 -l)(t-to-T). 

The numerical calculations confirm this estimation. The 
temperature T c and density p c of cold phase are determined 
by the condition of thermal equilibrium H = p c A(T c ) and 
the isobaric relation with hot phase, which corresponds to 
the relation p c = (p) [(1 - rj)T h /T c + rj\. 

The final mass fraction of the phases is an important 
characteristic of the models. We studied the dependence of 
the mass fraction 77 on the parameters A and to/t a . Fig. 4 
shows 77 as a function of to/tc for three values of the initial 
amplitude A = 10~ 2 , 10~ 3 10" 4 . In the case of a slowly ex- 
panding limit to /t c <; 10 4 the mass fraction 77 is a slowly in- 
creasing function of to/t c and weakly (< 5%) depends on A, 
when 77 iJ 0.85. In case of a fast expansion to/tc S5> 10 3 , the 
mass fraction 77 decreases, and its dependence on A becomes 
significant. The value of 77 reaches the minimum « 0.4 at 
to/tc w 10 2 , when the end time of the linear stage r ~ 0.3to 
[see equation (4.5)] becomes equal to the expansion time 
V/V ~ to/3 [see equation (4.1)]. For a faster expanding 
medium, to/t c Ss 10 2 , the final two-phase structure is not 
developed. In this case the expansion factor suppresses the 
thermal instability, and the models show the initial growth 
of the dispersion erf- until some maximum and the following 
decrease. The medium remains hot and uniform. 




I 1.02 1.04 



t/t 

Figure 5. The evolution of relative perturbations of tempera- 
ture or and pressure op for Hl-model and H2-model [solid lines, 
labeled as (1) and (2), respectively]. The evolution of temperature 
perturbations for the S-model (as in Fig. 2) is shown by a dashed 
line for comparison. 

The reconstructed ID space distributions of density and 
perturbed velocity for one period < x < A are shown 
in Figs 6 and 7 as dashed lines. The lines are presented 
at subsequent moments, and the period A increases with 
time as A(t) = Xot/t . The used reconstruction procedure is 
discussed in Section 3. 



5.2 Hydrodynamical models 

To solve equations (2.1)-(2.4) in the one-dimensional case we 
used an implicit hydrodynamical code. The code is based on 
the Lagrangian fully conservative numerical scheme of the 
second order accuracy (Samarskij & Popov 1980). 

We discuss here two models, which differ in the ini- 
tial wave length Ao of perturbations, Ao = Ap (Hl-model) 
and Ao = 5Ap (H2- model). They have the initial parame- 
ters A = 10 -4 and to/t c = 10 4 , which are identical to the 
S-model (see previous Section 5.1). In addition, we assume 
that the initial perturbed velocities and pressure variations 
equal to zero. The computed characteristics of the S-model 
almost coincide with those of the Hl-model. However the 
H2-model shows significant differences from both S and HI 
models in the linear and nonlinear stages. Fig. 5 shows the 
relative perturbations (square root of the correspondent dis- 
persion) of temperature <tt and pressure up as functions of 
time for HI and H2 models (solid lines). In the same Fig. 5 
the evolution of <jt for the S-model is shown for comparison 
(dashed line). The differences between the curves (Tr(t) for 
HI and H2 models are explained by the time delay. This time 
delay arises in the linear stage of evolution due to differences 
of the growth rates of different modes (see Appendix B, Ta- 
ble Bl). The perturbations for the shorter wave models grow 
more rapidly, and the most rapid growth rate is observed for 
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the isobaric S-model, which corresponds to the limit A — > 0. 
The shape of <7t curves for HI and H2 models are very sim- 
ilar in the nonlinear stage (<jt > 0.1), but delayed in time. 

To study the problem of nonlinear interaction of the 
modes with different Ao we calculated the number of models 
where the initial conditions were taken as a superposition 
of two modes, Ao = Ap and Ao = 5Ap. In the case of equal 
initial amplitudes of the modes the long wave mode was to- 
tally suppressed by the short wave mode in the nonlinear 
phase. The long wave mode suppresses the development of 
the short wave mode when the ratio of the initial amplitudes 
<; 50. This result quantitatively well agrees with the estima- 
tions made in the linear theory (Appendix B). Indeed, from 
Table Bl one can find the difference of the exponential fac- 
tors of these two modes AN = 9.56 — 6.0 = 3.56. It gives us 
the ratio of amplitudes exp(3.56) = 35. 

The small value of relative pressure perturbation ap 
10 -2 ) is an indicator of the isobarity of the process. Both 
HI and H2 models have peak-like up{i) curves (see Fig. 5), 
and the peaks are located at the end of the linear evolution 
stage, when the most rapid growth of temperature perturba- 
tions takes place. In the case of the Hl-model the maximum 
of ap has a rather small value <Jp max — 0.005 and even 
at the nonlinear stage the process is almost isobarical. This 
value of apmax can be compared with the correspondent 
value followed from the linear theory, ap max ~ a — 0.006 
(see Appendix B, Table Bl). After the peak ap decreases 
rapidly. At this stage the sonic waves are generated during 
the contraction of the cold phase. These waves produce os- 
cillations in ap clearly seen in Fig. 5. The intensity of the os- 
cillations also rapidly decrease with time. In the case of the 
H2-model the relative perturbations of pressure reaches a 
rather high value ap max ~ 0.1 (the linear theory gives 0.067, 
see Table Bl) and slowly decreases in oscillation regime after 
the maximum. 

The density and velocity distributions for the Hl-model 
are presented in Figs 6 and 7 (solid lines) at subsequent mo- 
ments. In the same figures the correspondent distributions 
for the isobaric S-model (dashed lines) are also shown. One 
can see a very good coincidence between these two models. 
In Figs 6 and 7 the curves for the Hl-model are used with 
the time shift 8t c with respect to the moments for the curves 
for the S-model. This shift follows from the time delay aris- 
ing during the linear stage of evolution of both models (see 
Appendix B). The shift is clearly seen as a small divergence 
of the <tt curves for the S and HI models in Fig. 5. 

The density distributions of the H2-model are shown 
in Fig. 8 at subsequent moments. The qualitative difference 
of density distributions in the HI and H2 models consists of 
the origin of 'wave hunches' in the latter case, when the cold 
phase is compressed by the hot phase. The wave hunches 
have a time period of about 27rf c (Ao/Ap). The origin of 
wave hunches indicates that the compression zone of the 
cold phase at the nonlinear stage is restricted by the zone 
of sonic wave propagation. 

The final mass fractions of the cold phase r] = 0.86 and 
0.885 for the hydrodynamical HI and H2-model are shown 
in Fig. 4 by a pentagon and square, respectively. The longer 
wave H2-model has the larger value of r\. 



10 - 
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x/A 

Figure 6. The density distributions at subsequent moments 
t/to = 1.0290, 1.0300, 1.0305, 1.0311, 1.0321, 1.0339 and 1.0560 
for the isobaric S-model (dashed lines) and hydrodynamical Hl- 
model (solid lines). The density distributions of these models are 
almost identical through the evolution. The full circles show the 
locations of the mass coordinate, which will finally divide cold 
and hot phases. 
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Figure 7. The distributions of velocity at subsequent moments 
(labeled by the numbers) t/t = 1.0290, 1.0300, 1.0305, 1.0311, 
1.0321 and 1.0339 for the isobaric S-model (dashed lines) and 
hydrodynamical Hl-model (solid lines). 
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Figure 8. As in Fig. 6, but for hydrodynamical H2-model. The 
origin of 'wave hunches' in the cold and dense phase, which are 
not seen in Fig. 6 for the Hl-model, indicates the breaking of the 
isobaric condition. 

6 DISCUSSION AND CONCLUSION 

We have developed a new statistical approach to study the 
thermal instability in an optically thin plasma. We inves- 
tigate the time evolution of the mass distribution function 
over temperature f{T), which characterizes the statistical 
properties of the medium. This medium can contain arbi- 
trary spaced three-dimensional perturbations. We construct 
the time evolution equation of function <p(T), which de- 
scribes perturbations of arbitrary amplitude and is correct 
under the isobaric condition. Formally, the obtained evolu- 
tion equation does not explicitly depend on the length scale 
A of perturbations, but the isobaric condition restricts the 
characteristic length of perturbations, A Ss Xp = 2nc 3 t c - 

The reconstruction of the space distributions of the den- 
sity, temperature and velocity field from the given tp(T) is an 
incorrect procedure (but we can recover these distributions 
under the restricted assumptions of the temperature profile). 
However, the statistical approach may be useful to give us 
the average parameters of the inhomogeneous medium, such 
as the pressure P = j{p)T, the total L = J p 2 A(T)dV = 

( P )MTT±) and spectral L e = J p 2 K{T)dV = (p)MT(^) 
luminosities. 

We study the evolution of function <p(T) numerically 
in the frame of the slowly expanding cloud scenario (see 
Section 4). We follow the thermal instability process in the 
linear and nonlinear stages including the formation of a 
two-phase medium. In the calculations we use the cooling- 
heating function (5.1), which does not allow the existence of 
a two-phases equilibrium. We find that in the case of the ex- 
panding medium the mass fractions of the phases have final 
nonzero values, contrary to the case of the stationary (not 
expanded globally) medium, where all mass asymptotically 
tends to the cold phase. We also find that the development 
of a two-phase medium can be suppressed by a fast expan- 



sion, when to/tc < 10 2 . The final mass fraction of the cold 
phase is larger for the slower expanded medium. Also, the 
final mass fraction depends on the amplitude of initial per- 
turbations. Usually, the larger amplitude produces a larger 
final mass fraction of the cold phase. These results of the 
calculations are general with respect to the assumed geom- 
etry of perturbations. In the case of one-dimensional and 
periodical perturbations we compare them with the results 
of hydrodynamical calculations. We show that our statistical 
isobaric approach correctly describes the thermal instability 
process on the length scales up to A ~ Xp. The hydrody- 
namical models show a weak dependence on the final mass 
fraction on the length scale in case of A > Xp. 

The presence of a quite strongly tangled magnetic field 
of character length scale A m <C Xp can suppress the thermal 
instability, because the magnetic pressure will obstruct to 
the contraction of the cold phase. In the case of a large scale 
magnetic field A m ^ Xp the contraction of cold gas will occur 
mainly along magnetic strength lines forming thin and flat 
structures (see also Meerson & Sasorov 1987 and Sasorov 
1988, where the formation of flattened plasma condensations 
was studied). 

In the case of a large scale magnetic field A m <; Xp 
the influence of the electron thermal conduction can be im- 
portant. Our approach is valid in the limited time interval, 
to < t < to + t + t K , when the hot phase temperature is 
not too high, Th < T K . At temperature Th — T K the heating 
rate (5.2) of the hot phase is compensated by the cooling 
effect of the thermal conductivity Th = K,T K p/ RphX 2 , where 
k = aT 5 / 2 is the coefficient of thermal conductivity (Spitzer 
1962), and a = 10 -6 in cgs units. In case of t K <C to the esti- 
mation gives T K ~3-10 5 {X/X P ) 4/7 K, and t K ~ (T K /T )t c ~ 
30(A/A P ) 4/7 t c . At short scales A <, W~ 2 X P the thermal con- 
duction suppresses the thermal instability even in the linear 
stage (Field 1965). The perturbations of length A ~ Xp can 
grow up to the maximum temperature ~ 3 ■ 10 5 -K\ Note 
that these estimations do not depend on the density and 
heating rate. The thermal conductivity facilitates the mass 
exchange between cold and hot phases. The following evo- 
lution of the two phases medium will be determined by this 
mass exchange process (see Aranson et al. 1993; Aharonson 
et al. 1994). 

There is not a full understanding of the problem of dom- 
inate scales of perturbations resulting in the thermal insta- 
bility process in the smooth scenario. It is known that the 
fastest growth rate has the modes with A & Xp. But, as it is 
shown earlier, the electron thermal conduction can suppress 
the growth of the shortest wavelength modes. On the other 
hand, the longer wavelength modes, A > Xp, have a slower 
development with an increase of the wavelength. Thus, there 
must be a dominate scale, which is determined by the com- 
petition of these factors. Another important factor, which 
has an influence on the value of the dominate scale, is the 
initial (at the moment to of the marginal stability) spec- 
trum of perturbations. This spectrum depends on the past 
history of the medium. We can just note that the linear the- 
ory predicts the dumping of any perturbations at this stage, 
t < to, and that the shorter wavelength modes are dumped 
faster. Finally, it seems that the dominate scale will be lo- 
cated close to Xp for a typical condition of the optically thin 
astrophysical plasma. 

The possibility of the two-phase medium model for 
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quasars was first discussed by McCray (1979). The proposed 
expanding cloud scenario of thermal instability can be rele- 
vant to the formation of this two-phase medium, when warm 
(T ~ 10 4 K) clouds are formed as the result of thermal 
instability in an expanded gas fragment, which is a prod- 
uct of either a star-star collision (Spitzer & Saslaw 1966) or 
star-accretion disc collision in the ~ lpc region around the 
quasar. The gas would be heated and ionized by the UV and 
X-ray emission of the quasar. 




(7-1) 
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APPENDIX A: NUMERICAL METHOD 

We briefly describe the numerical procedure to solve the 
equations (3.12) and (3.15) of the statistical approach. The 
equations are solved on the interval (Tmin, Tmax) of definition 
of the distribution function tp(T). The interval is divided on 
iV subintervals by the grid points T, where i = 0, 1, • • • , N, 
and To = T m i n , Tn = T max . The evolution of each point 
Ti is defined by the equation (3.12). If the distribution of 
{T™} are known at moment t n , one can find the new values 
of {T™ +1 } at moment t n+1 solving the system of algebraical 
equations 



At™ 
where 



T" 1 

— - 2 \ F i 



+ F? 



i = 0,l, 



(Al) 



(?) 



JV-l 

E 

i=0 



ACT) 
T 



dT. 



Here At" = t n+1 - t n is the time step. The upper in- 
dices denote the corresponding time. The function (p k (T) is 
piecewise-linear distributed through its average grid values 
T k ), where 



fi+l/2 - *i+l/2 



/Phi 



-1/2 



(p(T)dT, 



1 = 0,1, 



(A2) 



are conserved quantities, $ i+1 /2(i) = const, as it follows 
from equation (3.15). The algebraical equations (Al) are 
solved by the Newton-Raphson method. The time step At™ 
and the number of iterations are determined by the assumed 
accuracy and convergency of the solution. Described numer- 
ical procedure can not be applied to calculate the evolution 
stage of highly inhomogeneous medium, when the time step 
becomes very small due to small differences Ti+i — T. In 
this case we modify the numerical algorithm introducing a 
grid redistribution procedure at each time step. Namely, we 
choose new values of {T} inside of the interval (T m i n , T max ) 
and remap {^i+1/2} from the old values to the new ones. 



APPENDIX B: LINEAR THEORY IN THE 
SMOOTH SCENARIO 

The linear theory of thermal instability (see Weymann, 1960, 
Field, 1965) provides us the relations 

Pi 



_ _ Q Pi 

Po po 



a Ti 
1 + aT? 



(Bl) 



y£ x , between the perturbed 



where a = ■y(\n/2irc s ) 
pressure Pi, density p± and temperature T\ with a wave- 
length A = £Xp, also the velocity variation v/c 3 = ixlpi/po, 
and also the growth rate n — x jt c of the condensation mode, 
where a; is a positive real root of the equation 



e 2 x 3 + , 



i 2 x 2 + X - 



x . 



(B2) 



Other two roots correspond to the dumped wave modes. 
Here xq = not c , /3 — 7(7 — 1). The scale parameters Ap, 
t c and no are given by equations (3.1), (3.2) and (3.17'), 
respectively. We supply here and below the index to all 
relevant parameters (such as n, N, r), which are used in the 
main text [see formulae (3.17'), (4.3)-(4.5)]. 

In a frame work of the smooth scenario (Section 4), 
when £ = tc/to -C 1, the parameter xo(t) increases slowly 
and changes the sign from negative to positive at t = to 
linearly [equation (4.3)] 

, . St 
xo(t) = a — , at — t — to, 
to 



(B3) 
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where a = 3b(7 — 1). In this scenario the linear theory at 
St > gives the following: 

1. The £ = case corresponds to the isobaric S-model. We 
set I = into equations (Bl) and (B2), and find a = and 
the root x — xo. The linear stage of the S-model is described 
by the equations (4.3)-(4.5) for the no = Xo/t c , the integral 
No(t), and the end time to, respectively. The integral 

xl{t) 



N (t) 



2< 



(B4) 



reaches the No — Lq — \n(l/Ao) level (when the amplitude 
A = Aoe L ° = 1) at the end time 

n=X - = -y/2aLoe, (B5) 
a a 

when the end (maximal) growth rate reaches 

Xo = xo (t ) = A/2aLoe (B6) 

Note that x < X ~ £ 1/2 < 1, when £ < 1. 
2. The I > cases correspond to the hydrodynamical models 
(i = 1 for the HI and £ = 5 for the H2-models). The integral 
Ni has an implicit form now 

V, = - / xdt=^(l + -pfx+-£ 2 x 2 ), (B7) 



tc 



where we use dt — todxo/a [see equation (B3)], dxo/dx fol- 
lows from equation (B2) and x is defined by equation (B2). 
The end rate Xe can be found from equation 

2a£Le = Xf(l + ^f3£ 2 X e + (B7') 

which is followed from equation (B7) where we set the end 
value Ne = Lit = \n(l/A oe ). Note, that Xe < sflaLdl ~ 
^1/2 ^ -y -pjjg enc j time r£ follows from equations (B2) and 
(B3), 

n = -{X t + P£ 2 Xj + £ 2 Xf), (B8) 
a 

where Xt is taken from (B7'). 

3. It is useful to compare No and Ne [equations (B4) and 
(B7)], because the difference AN = No(t)—Ne(t) determines 
the suppression of the hydrodynamics modes {I ^ 0) relative 
to the isobaric mode (I = 0). Substituting xo from equation 
(B2) to (B4) we find the difference 



where x is determined by (B2) and 
3a; , 



J + ^[l + 2£ 2 (/3 + x) 2 ]. 



(B9) 



(BIO) 



It is also useful to compare the end times re and fo = 
To(LefLo) 1 ' 2 [see equation (B5)], where to is the end time 
of I — mode when the initial amplitude of this mode is 
equal to Aoe- The end times difference can be obtained from 
the form 



r?-f° 2 = l(^yW 2 Xtg(X e ), 



(Bll) 



where Xt is taken from equation (B7'). 
4. Note that in our case ( < 1, when any x, xo, Xo, Xe 
~ ^Z 2 <^ 1, the terms x 2 in the brackets in equations (B7) 
and (B7'), the terms x 3 in equations (B2) and (B8), and the 



Table Bl. Parameters of the S-model (I = 0), Hl-model (£ = 1) 
and H2-model (£ = 5) followed from the linear theory. 



i 


x(t ) 


N e (ro) 


x e 


Tl /to 


a(X e ) 





0.063 


10 


0.063 


0.0316 





1 


0.059 


9.56 


0.060 


0.0323 


0.006 


5 


0.033 


6.0 


0.040 


0.042 


0.067 



i-term in the bracket in equation (B10) can be omitted. The 
equation (B2) is a quadratic now: f3£ x + x = xo (t) and has 
a root 



= x(t) = (jl+4:/3£ 2 Xo(t)-lj /2(3£ 2 



(B12) 



Substituting x(t) in equations (B7) and (B9) we obtain the 
explicit time-dependencies of Ni(t) and AN(t). We also find 
a(t) [see equation (Bl)] and can now construct both the 
a T (t) = a e N(t) and a P (t) = a(t)a T (t)/(l + a (t)) functions 
on the linear stage, at < St < T(. Graphs of these functions 
coincide well with the left-hand branches (at moments < 
St < rt) of the corresponding lines in Fig. 5. 
5. In the case £ 2 xo <S 1 the above written formulae can be 
reduced to the simpler forms. From equations (B2, B12) we 
find 

x ~ x - (3£ 2 x 2 o, x < xo <S -4-. (B13) 
The difference (B9) in the main (x = xo) order is 



AN(t) ~ ^xg(t). 



(B14) 



This difference at the end time to of the £ = mode, when 
xo = Xo ~ \/2aLo£, [see equation (B6)], is 



ATV(to) ~ -/3fLo^2aL t 

The end rate Xe can be found from equation (B7') 

X e ~ y/2aL£ (l - ^f3£ 2 y / 2aL e ^ . 

The end time [equation (B8)] is now 

n ~ ^ ^2aUi (l + \pf^2aL t £\ , 

and the end times difference [equation (Bll)] is 

ti> — to~ ^-(3£ 2 Let c . 



(B14') 



(B15) 



(B16) 



(B17) 



6. For numerical examples we take b = 1, 7 = 5/3, a = 2, 
(5 = 10/9, £ = t c /t = 10" 4 , L = he = 10. Results are 
shown in Table Bl and can be compared with the results of 
the statistical theory and exact hydrodynamical calculations 
presented in Section 5. 

This paper has been produced using the Royal Astronomical 
Society/Blackwell Science style file. 
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